Dynamics and hysteresis in square lattice artificial spin-ice 
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Dynamical effects are considered in a model for a spin ice on a square lattice in two dimensions. 
The spin ice is an array of thin elongated magnetic islands interacting with each other via long-range 
dipolar forces, which leads to geometrical frustration. Each island is assumed to possess a three- 
component dipole moment free to point in any direction, but with strong uniaxial anisotropy energies 
that select preferred directions. The model has real dynamics, including rotation of the magnetic 
degrees of freedom, going beyond the Ising-type models of spin ice. The dynamics is studied based 
on a Langevin equation that includes the temperature and a phenomenological damping, which is 
solved using a second order Heun algorithm. As an important application, the hysteresis in an 
applied magnetic field is calculated for typical island sizes used in recent experiments. At the same 
time, thermodynamic properties such as the specific heat can be studied - a peak in specific heat is 
related to a type of melting-like phase transition present in the model. 
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INTRODUCTION: SQUARE SPIN ICE, 
FRUSTRATION, DYNAMICS 



Artificial spin ices are systems in two dimensions that 
mimic the usual three-dimensional spin ice materials that 
exhibit geometrical frustration effects: not all the pair- 
wise spin interactions can be satisfied simultaneously^—. 
Indeed, the name spin ice comes from the fact that these 
systems obey the ice rule: for a square lattice, in each 
vertex four spins meet but two must point inward while 
the other two must point outward. In general, the ar- 
tificial compounds are built from magnetic nanoislands 
(generally made of permalloy) which can be organized in 
several different geometries in such a way that the frus- 
tration is manifested^—. 

Here, the focus is on the artificial square lattice spin ice 
array, first fashioned and studied by Wang et al£ in 2006. 
This artificial square ice consists of magnetic nanois- 
lands (with a shape that looks like a "cigar" ) arranged as 
shown, for instance, in Ref. H|. Each nanoisland contains 
a net magnetic moment that tends to point out along the 
longest axis. When the interactions between neighboring 
islands are increased, the system of Ref. was increas- 
ingly filled with vertices that obeyed the ice rule (two- 
in/two-out, which is non-degenerate in two dimensions). 
Despite this, the predicted ground state of square ice was 
not observed experimentally until the work by Morgan et 
alM- in 2011. Using more advanced techniques, those au- 
thors could observe large regions of their samples which 
were able to adopt the square ice ground state. With 
this progress, they could also observe the predicted exci- 
tations above the ground state, which resemble magnetic 
monopoles connected by energetic strings^— (similar 
to Nambu monopoles^ 7 -). These elementary excitations 
are different from those of natural three-dimensional spin 



ices, which are magnetic monopoles connected by observ- 
able but non-energetic strings^. Therefore, these arti- 
ficial compounds have attracted large interest in recent 
years. 

The manipulation and control of the excitations could 
lead to several technological applications. However, there 
is a primary problem with these artificial systems, par- 
ticularly with the square ice: the lack of thermaliza- 
tion. Thermal activation of the island's magnetic mo- 
ment ("spin") configurations is very weak or nonexistent. 
Lack of thermalization is an important topic for the ex- 
perimental artificial spin ices and it can be partially al- 
leviated by applying varying external magnetic fields 7 -^. 
Moreover, reductions in island volume and magnetic mo- 
ment through state-of-the-art nanofabrication can bring 
energy scales closer to room temperature, leading to ther- 
mally driven slow dynamics. Alternatively, the use of ma- 
terials with an ordering temperature near room temper- 
ature seems to be another important possibility. Really, 
by using such a material, a recent experimental work on a 
square lattice in an external magnetic field confirms a dy- 
namical "pre-melting" of the artificial spin ice structure 
at a temperature well below the intrinsic ordering tem- 
perature of the island material, creating a spin ice array 
that has real thermal dynamics of the artificial spins over 
an extended temperature ranged. Another expectation 
for better understanding of these compounds is to study 
similar systems such as colloidal systems, which have an 
advantage over the usual magnetic arrays because ther- 
mal activation of the effective spin degrees of freedom 
is possible^. So, a more detailed analysis of the effects 
of thermal fluctuations and the spin dynamics in a two- 
dimensional spin ice material should be of great interest 
for a better understanding of these interesting frustrated 
systems. 
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Particularly, thermal effects in artificial square ice 
were studied recently by some of us^S. The focus was 
to examine the roles of elementary excitations in the 
thermodynamic properties of these systems. We have 
found that the specific heat and average separation be- 
tween monopoles with opposite charges exhibit a sharp 
peak and a local maximum, respectively, at the same 
temperature^, T p w 7.2D/ks, where D is the strength of 
the dipolar interaction and ks is the Boltzmann constant. 
These features increase logarithmically with the lattice 
size and, therefore, the possibility of string breaking and 
free monopoles were discussed. However, all calculations 
were done by assuming an Ising-like behavior for the net 
magnetic moment of the nanoislands. The Ising behav- 
ior of the islands seems to be the most realistic situation 
for the typical artificial magnetic ices made of permalloy 
and therefore, these systems could not display real dy- 
namics; it is one of the causes concerning the problem of 
lack of thermalization. In this article, we study possibil- 
ities beyond the Ising behavior for the nanoislands and 
consequently, our attention shifts to magnetic ices with 
real dynamics and the extra features that such dynam- 
ics may produce. For this study, the internal structure 
and shape of the magnetic nanoislands have to be taken 
into account in details. The theoretical study of the net 
magnetic moment (with more degrees of freedom) of in- 
dividual magnetic nanoislands with different sizes and 
shapes is the initial point. Such an investigation could 
also predict better circumstances for novel ideas and dic- 
tate protocols for the realization of artificial magnetic 
spin ices with real thermal activation and spin dynamics. 

A detailed study concerning non-Ising spins for indi- 
vidual islands can be found in Refill. By applying the 
techniques and ideas of our previous paper—, we now 
consider a more intricate situation concerning artificial 
square spin ices. Here we study the dynamics of the 
array by assuming that a nanoisland's spin is free to 
point in any possible direction, but with strong uniaxial 
anisotropy energies that favor preferred directions. Then, 
differently from all previous articles published until now 
on this topic, which consider basically only the dipolar in- 
teractions among the islands, here, additional terms con- 
cerning these anisotropics must be included in the Hamil- 
tonian. We would like to answer questions like these: 1) 
Is there any possibility of building artificial magnetic ices 
(for instance, by using permalloy islands) exhibiting real 
dynamics, only by decreasing the island sizes? 2) If there 
is this chance, what are better parameters for the islands 
and for the arrays with spin dynamics? 3) If there is not 
such a chance, what should the parameters be for other 
potential magnetic materials and their ice arrays, so they 
exhibit real dynamics? 

The article is organized as follows: In section II, the 
model is explained in detail. We also define two order pa- 
rameters, which are important to identify if the ground 
state of the models analyzed can be accessed at low tem- 
perature. Five models (denoted by A, B, C, D and E) 
with different lattice and island parameters are proposed 



to see the possibility of a spin ice dynamics. The nu- 
merical calculations are presented in Section 777; The 
dynamics is investigated by using a Langevin approach 
for the spins. In Section IV, some thermal equilibrium 
properties of the models A, B, C, D and E are calculated. 
In Section V, we present the hysteresis calculations and, 
finally, the conclusions and some discussions are given in 
Section VI. 



II. THE MODEL SYSTEM 

The open square ice system with N c = L\ x L2 unit 
cells can be set up as follows. One can define the sites 
of a square lattice in the usual way: the k th lattice site 
(a monopole charge center) is at a point r k = (xk,Uk), 
where x k — rrika and y k — n k a are integer multiples of 
the chosen lattice constant a. The points are chosen to fit 
inside the desired L\ x L2 area. For each unit cell of size 
a x a there are two nano-islands that act as a two-atom 
basis, with the locations, 



rki = (m k + \ ,n k )a, 
rk2 = (mk, n k + \)a, 



(1) 



where the "1" and "2" refer to the sublattices. The nano- 
island on the 1 st sublatticc has its long axis in the x 
direction; the other nano-island, on the 2 nd sublattice, 
has its long axis in the y direction. There are N = 2N C 
islands in the whole system. 

A 3D vector magnetic moment jli, i = 1,2, 3...N is as- 
sociated with each island, whose defined center position 
is some fi . Each f*j is selected from the set of ffc CT , where 
(7 = 1,2 denotes the sublattice. We use indeces k, I for 
locations of the unit cells, and indeces i,j for locations 
of individual islands or their dipoles. 

The dipoles are assumed to have fixed magnitude /x, 
while their direction is represented by a unit vector jli. 
The magnetic moments interact via long-range dipole 
forces, and are also affected by two forms of local shape 
anisotropy. First, there is a uniaxial anisotropy that im- 
pedes free rotation in the xy-plane, associated with some 
energy constant K±, and oriented along x for the first 
sublattice and along y for the second. Dependent on its 
sublattice, each moment has an axis Ui (equal to x or y) 
for this anisotropy. Second, because the nano-islands are 
thin in the z-direction, the z-direction is a hard axis, and 
there is a hard-axis anisotropy whose energy scale is de- 
termined by a constant K%, the same for all the islands. 
The Hamiltonian is then 



^ = fig [I 2 [3(jlj ■ fjjXAj ■ fjj 



(2) 
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i 

Here /iq is the magnetic permeability of space, and fij is 
the unit vector pointing from the position of jlj towards 
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the position of fli. The first sum is the dipole-dipole 
interactions, the second sum contains the anisotropy en- 
ergies and an applied external magnetic induction B ext = 
/io-ffcxt- A constant is included in the K% anisotropy en- 
ergy so that that energy is zero when a dipole points 
along its local anisotropy axis u$. It is interesting to note 
that if a dipole moves in the xy plane, it only pays the 
cost of the K\ anisotropy term, but motion up out of 
the xy plane (say, in the xz plane) involves an energy 
proportional to the sum of both anisotropics, K\ + Kj,. 

The motion out of the xy plane is also impeded by 
the dipolar interactions. With the dipole pair distances 
scaled by the lattice constant, the effective strength of 
nearest neighbor dipolar interactions is determined by 
the dipole energy factor, 



47r a 3 



(3) 



Depending on the island geometry, which is discussed 
further below, the anisotropy constants K\ and K3 
would typically be of similar order of magnitude. Thus, 
there are three important energy scales: dipolar energy, 
anisotropy energy, and the thermal energy ksT. The 
anisotropy constants are proportional to the volume of 
the islands, as is /1 = M S V , where M s is the saturation 
magnetization of the magnetic material. But then, this 
dipolar constant D increases as the squared island vol- 
ume. Thus, it is possible to adjust the island size (and 
spacing a) to get these energy scales as desired in relation 
to each other. Typically, the interesting case must have 
the thermal energy less than both the effective dipolar 
energy (per site) and anisotropy energy. But note, the 
effective dipolar energy can be quite a lot larger than 
that indicated by D, which only measures the energy in 
a nearest neighbor pair. When the dipolar interactions 
are summed, the net dipolar energy associated with any 
island could be much larger than D. 



A. Spin-ice ground state and order parameters 

For the square lattice spin-ice, the ground state is two- 
fold degenerate, and involves alternating dipoles on each 
of the two sublattices. The ground state fully satisfies the 
two-in/two-out rule in each monopole charge cell (junc- 
tion of four islands at the site ffc of each unit cell) . The 
unit cell positions are expressed = (rrik,nk)a, where a 
is the lattice constant and rrik and rik are integers. Then 
one of the ground states can be constructed by setting 
the dipole directions as: 



<:s - Af s H) = +(-ir fe+nfc * 



A?l = $ S (r k ) = -(-1)™*+"* 



x, 



(4) 



This formula is arranged so that at a chosen unit cell at 
position ?k, the dipole on one sublattice points inward 
and the dipole on the other sublattice points outward, 
thereby globally enforcing the two-in/two-out rule. By 
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FIG. 1: A 16 x 16 model system with d = fci = fc 3 = 0.1, 
in a metastable state at temperature fcsT/e = 0.025, from a 
hysteresis scan (this is a state at h ex t = 0). Most of the system 
is locally close to the Z = +1 ground state. The upper right 
hand corner is locally near the Z = — 1 ground state, and there 
is a bent domain wall connecting the two regions. For interior 
charge sites (junction points of four islands), there happens 
to be no discrete monopole charges present: all — and 
the discrete p m = 0. 



reversing the sign on all the dipoles, the other ground 
state is obtained. 

With the ground state determined, we can construct 
a measure of the proximity (in phase space) of any arbi- 
trary state to one of the ground states. This order pa- 
rameter Z is simply the overlap with this ground state: 
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(5) 



k=l ir=l 



The index a labels the sublattice. If the system happens 
to be found in the ground state defined in Eq. (j4]), then 
Z = 1; if the system is in the inverted ground state, then 
Z = — 1. Thus it is possible to show that the range of 
Z is from -1 to +1. This order parameter is useful for 
indicating the degree of thermodynamic excitation in the 
system, by the deviation of \Z\ from unity. Further, its 
sign then gives an indication of processes which involve 
the transformation from one ground state to the other. 
Indeed, considered even as a local variable (calculating 
only near a single charge cell), we can track when the 
system has different regions close to either of the ground 
states, possibly with regions separated by domain walls. 
Fig. [1] shows an example, where the system has a region 
near one of the ground states, with Z « +1, separated by 
a domain wall from another region that is near the other 
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ground state, with Z f» — 1. The net averaged value of Z 
for the entire system, however, acquires an intermediate 
value, Z w 0.521, indicating considerable separation from 
a uniform ground state. 

The other obvious order parameter to be measured is 
the areal density of monopole charges, p m . We make a 
simple discrete definition, to connect to Ising spin ice 
models, and, a more generalized continuous definition 
that accounts for the greater freedom of the continuous 
dipoles in the model described here. The discrete defi- 
nition of a monompole charge involves counting the net 
number of dipoles that point outward at a chosen charge 
site ffe, and dividing that result by two. There are four 
dipoles fii k ,ik = 1,2,3,4, surrounding any charge cell 
center r^. Then the possible monopole charge values are 
qk = 0, ±1,±2, although the last values may typically 
be of low probability. For the discrete charge definition, 
whether a dipole points outward or inward is determined 
with a Heaviside step function H(x): 



qk 



i 4 

? £[2#(/V«iJ-l] 



(6) 



The unit vectors ii k ,ik — 1,2,3,4, point outward from 
charge site r^ to each of the four nearest islands. 

Because this discrete definition can show sudden 
change when a dipole rotates 90° from the radially out- 
ward direction, we also considered a continuous defini- 
tion. In the continuous definition, the step function is 
removed, and only a scalar product is needed, 



Qk 



1 v- 



(7) 



ifc=l 



In contrast to the discrete definition, this charge defini- 
tion varies continuously from qk = — 2 to qk = +2. One 
can also note that for either the discrete or the continuous 
definition, total monopole charge is conserved. A posi- 
tive contribution produced by some dipole at one charge 
cell is accompanied by an equal negative contribution at 
a neighboring charge cell (each dipole contributes to two 
charge cells). Then, the total algebraic monopole charge 
in the system takes the conserved value, zero. 

In order to get a measure of the monopole charges 
present, regardless of their sign, we define a density for 
the system as a whole, by applying absolute value. Thus, 
the "monopole density" measured in the simulations here 
is defined as 



1 



k=l 



Qk\ 



(8) 



This is normalized by the number of charge cells. By 
using absolute value, the definition does not allow the 
cancellation of charges of opposite signs. Note that in 
either of the ground states, there are no charges at any 
sites, and p rn — 0. Charges appear as the system moves 



away from the ground state. Thus, this is another mea- 
sure of excitation in the system. (This is true for the spin 
ice on the square lattice, but not on the Kagome lattice, 
whose ground state contains charges, due to there being 
three dipoles for each charge cell.) 



B. The undamped dynamics 

The zero-temperature, undamped dynamics of each 
magnetic dipole is determined by a torque equation, 

dfli 
~dT 



(9) 



where Bi is the local magnetic induction acting on the i th 
dipole, and 7 is the electronic gyromagnetic ratio. The 
local magnetic induction is derived from the Hamiltonian 
by assuming an energy — fti ■ Bi for each dipole, i.e., 



Bi 



1SH 

M 6^ 



Ki. 
2 — (fa ■ Ui)ui 
A 1 



D 



3(/*j • Tij)Ti 



{rij/af 



2—(fi i -2t)S + B fat . (10) 
M 

It will be convenient to choose some standard units for 
the time, the applied field, and so on, to simplify and 
scale the numerical calculations. The dipole terms are 
simplified by selection of the lattice constant a as the 
unit of length. A natural unit to measure field H exi is 
the saturation magnetization M s from which the particles 
are made. For example, for Permalloy, with M s = 860 
kA/m, this unit as a magnetic induction is close to one 
tesla: /j,qM s « 1.08 T. Using this quantity to scale the 
magnetic field and hence the magnetic induction defines 
their dimensionless field, 



B„ 



HoM s 



(11) 



When h ex t approaches 1.0 the applied field should have a 
strong tendency to saturate the magnetization of the sys- 
tem (if the dipolar interactions do not impede that). This 
then indicates how to scale the dipole and anisotropy 
fields, i.e., by writing the dimensionless local magnetic 
fields from (fTDl). 



hi = 



Bi 



D 



E 



3(/ij ■ rij)r\j fij 



{rij/af 



(12) 
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K, 



-(Ai • z)z + h cxt . 



/j, Ho M s ppoMs 

This suggests the introduction of dimensionless coupling 
constants that indicate the relative strength of each con- 
tribution, 

D 



d = 



ft 



Ki . K 3 
, K3 = 



(13) 
(14) 
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These definitions involve the different energy scales di- 
vided by an energy unit, 

e = (J-o^Ms, (15) 

that depends on the size of the magnetic islands. The 
dimensionless dipole parameter d can be seen to be pro- 
portional to the volume fraction of the system occupied 
by magnetic islands, since /i = M S V for each island. Ob- 
viously a higher packing of magnetic material into the 
lattice leads to stronger dipolar effects, and d indicates 
their effective strength. 

The dynamic equation can be scaled in the same way, 
so that the dimenionless field appears on the RHS. Then 
the dynamics for the unit vector dipoles is described using 
a rescaled time r, 

^1 = x hi, t = 7/i Af s t. (16) 
dr 

With the above scaling of the fields, the unit of time 
is (7/ioM s ) _1 . For the case of Permalloy and using the 
gyromagnetic ratio as 7 = e/m e ps 1.76 x 10 11 T _1 s -1 , 
this unit is ('j^oM s )~ 1 ps 5.26 ps. 

C. Island geometry and energetics 

The shape anisotropy constants K\ and K3 can be es- 
timated based on the magnetic properties for Permalloy 
(or other material) and micromagnetics simulations for 
the choice of island geometries and island volume V. In 
this work we consider thin elliptical islands. Here L x 
denotes the major-diameter of the ellipse and L y is the 
minor-diameter, while L z is the height of the island or its 
thickness. The semi-major axis is A = L x /2, the semi- 
minor axis is B = L y /2. It is well-known that the an 
elliptically shaped magnetic particle will have a nonuni- 
form anisotropy 22 within the plane of the island. In Ref. 
I21I the anisotropy constants (as energies per unit volume) 
were estimated based on a calculational approach for thin 
elliptical islands, for a range of thicknesses L z <C L x , 
characterized by an aspect ratio #3 = L x j L z , and various 
lateral aspects ratio g\ = L x /L y . Here we consider some 
different sizes and shapes for the islands and discuss the 
advantages and disadvantages, as far as the expectations 
for their dynamics and and the relative importance of the 
different energy scales when placed in a square spin-ice 
array. 

Model A. In Wang et al&, experiments on square 
spin-ice were carried out for (quasi-rectangular) particles 
with dimensions 220 nm x 80 nm x 25 nm, where the last 
number is the vertical thickness. In those experiments, 
the particle sizes were kept fixed, but different lattice pa- 
rameters a from 320 nm to 880 nm were used. In this first 
model to consider, we use these numbers to describe ellip- 
tical particles: L x — 220 nm, L y = 80 nm, L z = 25 nm. 
Then the particle volume is V — nABL z = 3.46 x 10 5 
nm 3 , and using a saturation magnetization M s = 860 



kA/m for Py, the magnetic dipole moment per parti- 
cle is fi = 2.97 x 10~ 16 A-m 2 , the equivalent of about 
3.2 x 10 7 Bohr magneton. For its aspect ratio parameters 
gi — 2.75 and 173 = 8.8, the anisotropy energy densities 
can be found by interpolation of the simulation results in 
Refi2i, as K x /V = 0.0064 4/nm 2 and K 3 /V = 0.0143 
A ex /nm 2 , where A cx pa 13 pJ/m is the exchange stiff- 
ness for Py. The easy-axis anisotropy then works out to 
K 1 = 2.9 x 10~ 17 joules, while the hard-axis anisotropy is 
estimated as K 3 = 6.4x 10~ 17 J. These are considerably 
larger than room-temperature (300 K) thermal energy 
/c B T« 4.1 x 10~ 21 J, as needed for stable magnetic mo- 
ments. The energy unit is e = hqij,M s = 3.21 x 10~ 16 
J. Then the dimensionless anisotropies are k\ = K\/e = 
0.0897, k 3 = K 3 /e = 0.200 . The scaled thermal energy 
at room temperature is T = ksT/e = 1.29 x 10~ 5 , an 
extremely small value. 

The nearest neighbor dipolar energy scale might be es- 
timated first at lattice constant a — 880 nm, for which 
it is D = 1.29 x 10~ 20 J, about 2000 times smaller than 
K\. If instead the lattice constant a = 320 nm is used, 
this will scale up by a factor of (880/320) 3 , leading to 
D = 2.68 x 10" 19 J, or still 100 times smaller than K x . 
The dimensionless dipolar coupling for a — 320 nm is 
d = D/e = 8.35 x 10~ 4 . Obviously, values of ki, k 3 , and d 
similar to these are needed to get a spin-ice system, how- 
ever, dynamics simulations could be difficult with these 
parameters because the anisotropy is so dominant and 
Ising-like. Over the time scales that can be accessed 
in numerical simulations, one would not expect to see 
much dynamical flipping of the island dipoles, except in 
the presence of a strong applied magnetic field. Thus it 
may be interesting instead to consider some other parti- 
cle sizes where the dynamics can be expected to be more 
active. 

A thinner or smaller island will result in a smaller mag- 
netic dipole moment fi, which leads linearly to weaker 
anisotropy, but quadratically to weaker dipolar energy. 
Both energy scales become closer to the thermal energy. 
Thus we can try to change the particle size in such a way 
so that room temperature thermal energy is closer to Ki 
and perhaps even larger than D. 

Model B. Consider smaller elliptical particles, with 
L x = 100 nm, L y = 12.5 nm, L z — 5.0 nm, on a lattice 
with constant a = 200 nm. The aspect ratio parameters 
are gi = 8.0, g 3 = 20. In this case, the smaller particle 
volume V = 4910 nm 3 leads to a smaller magnetic dipole 
moment \i = 4.22 x 10 -18 A-m 2 , and the anisotropy en- 
ergy constants are found to be K\ — 6, 66 x 10~ 19 J and 
K 3 = 7.48 x 10~ 19 J. The energy unit is e = fiofiM s = 
4.56 x 10~ 18 J. This leads to scaled energy parameters 
fci = K x /e = 0.146, fc 3 = K 3 /e = 0.164, and for lattice 
constant a — 200 nm, one has d = D/e = 4.87 x 10~ 5 . 
The effective dipolar coupling is very small in this case. 
The scaled room temperature thermal energy (300 K) is 
now T = ksT/e — 0.908 x 10~ 3 , stronger than the dipo- 
lar coupling d. Thus this model might exhibit observable 
changes in its properties for easily accessible tempera- 
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tures. 

Model C. One can also consider even smaller par- 
ticles, with L x = 40 nm, L y = 8.0 nm, L z = 4.0 
nm, to try and get weaker anisotropy energy scales 
(still using the Py parameters). The particle volume 
is now only V — 1005 nm 3 , and the dipole moment is 
/i = 8.64 x 10~ 19 A-m 2 . At aspect ratio parameters 
gi = 5.0, #3 = 10, the anisotropy energies are found 
to be somewhat smaller: K\ = 1.38 x 10~ 19 J, and 
K3 = 1.12 x 10 -19 J. The energy unit, though, is now also 
smaller: e = 9.34 x 10 -19 J, leading to dimensionless cou- 
plings k\ — 0.148 and ^3 = 0.120, only slightly different 
from those in Model B. However, the smaller energy unit 
means that room temperature effects may be more acces- 
sible. The scaled thermal energy at 300 K is increased: 
T = k B T/e = 0.00443 . For a lattice with a = 80 nm, 
we find D = 1.46 x 10~ 22 J, and d = De = 1.56 x 10~ 4 . 

It is clear from the above examples, and also from the 
numerical results presented below, that the strong Ising- 
like anisotropy for real spin-ice particles dominates over 
thermal energy at (and below) room temperature. That 
being the case, we find it interesting also to study some 
models with fictitious parameters, which might be possi- 
ble to achieve in other materials with different values of 
M s , Kx/V, etc. 

Model D. Rather than assuming a particular particle 
size and using Py parameters, suppose some particles 
are arranged so that D = K\ = K3 = j^e. Obviously 
nature may not easily produce such a system with all 
equal energy scales, but it may be possible by appropriate 
materials engineering. We use a fraction of e, which is 
required by the definition of d, see Eq. (1131) (the volume 
fraction of dipoles on the lattice cannot be more than 
unity) . Then the scaled energy parameters are all equal: 
d = k\ = &3 = 0.1 . A physical value of e is needed, 
based on values of fi and M s for some real particles, to 
locate room temperature on the temperature scale. 

Model E. This is similar to Model D, but with a 
stronger value for the hard-axis anisotropy: D = K\ = 
jqS, and K3 = \e. This would account for the typi- 
cally stronger anisotropy for out-of-plane dipolar motion, 
due to the assumption of thin nano-islands. Only some 
limited simulations were done for this case, because the 
results are found to be very similar to those for Model D. 



III. NUMERICAL CALCULATION: LANGEVIN 
DYNAMICS 

The dynamics is investigated here using a Langevin 
approach for the "spins" fii. This requires including the 
combination of a damping term and a rapidly fluctuat- 
ing stochastic torque in the dynamics. The size of the 
stochastic torques is related to the temperature and the 
damping constant, such that the system tends towards 
thermal equilibrium for the chosen temperature. The 
approach also gives the dynamics in the limit of zero 



temperature, if needed: the stochastic torques are set to 
zero, but the damping can still be included. 

In fact, it is more physical to think that the dynam- 
ics includes random magnetic fields, rather than torques. 
The dynamical equation for some selected unit dipole ex- 
posed to a deterministic field h and a stochastic field h s 
is written in the dimensionless quantities as 



djj. 



fi x (h + h s ) — afi x 



(fi x ( h + h t 



(17) 



The first term is the free motin and the second term is the 
Landau-Gilbert damping, with dimensionless damping 
constant a. The dynamics evolves according to the su- 
perposition of the deterministic and random fields. One 
can also split out their contributions separately, 



dfi 

7h 



— = LI X 



h- 



(ji x hj 



+ fi x 



h s — a ( fi x h 



which shows the effective dynamical fields (in the brack- 
ets) acting on this dipole. Then in order for the stochastic 
fields to establish thermal equilibrium, they are assumed 
to have time correlations determined by the fluctuation- 
dissipation (FD) theorem, 



(K(r)hi(T'))=2aTS ij S(r-r'). 



(19) 



The indices i, j refer to any of the Cartesian coordinates. 
The dimensionless temperature T is the thermal energy 
scaled by the energy unit, 



T 



k B T 



k B T 



(20) 



The fluctuation-dissipation theorem indicates that the 
power in the thermal fluctuations is carried equivalently 
in the random magnetic fields. For reference, with the 
physical units replaced the relation can be expressed 

jfi(B i s {t)Bi(t'))=2ak B TS ij S(t-t'). (21) 

The Langevin equation in (|18p is a first-order differ- 
ential equation where the noise is multiplicative - the 
second term involves a product of fx with a net stochas- 
tic field. To discuss the solution method, it is simplest 
to let y = y(r) be a vector that represents the entire 
set of spins, y = {/^(r)}. Then symbolically y obeys a 
differential equation in the general form, 

^=f[r,y(r)] + fs[r,y(r)}-h s (r). (22) 

The function / represents the deterministic time deriva- 
tive on the RHS of (fT5|) and the function f s represents 
the stochastic part of the dynamics. Each is defined in- 
directly by comparing this with the Langevin equation. 
The fields /, f s and h s are vectors of 3iV components, 
where N is the number of dipoles in the array. 
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A. Second order Heun integrator 

An efficient method for integrating this magnetic dy- 
namics type of equation forward in time is the second or- 
der Heun metho d 23 ' 24 . That is in the family of predictor- 
corrector schemes and is rather stable. Another method 
that could be applied with similar precision would be 
second order Runge-Kutta, however, it would be slightly 
less efficient for including the stochastic fields. 

The predictor stage for the second order Heun algo- 
rithm is an Euler step, which is followed by a corrector 
stage that is equivalent to the trapezoid rule. Each in- 
volves moving forward in time over some time step At, 
with the needed results obtained by integrating Eq. [22l 
from an initial time t„ to a final time t„ + i = t„ + At, 
during which the stochastic fields are acting. With nota- 
tion y n = y(r n ), the predictor stage produces an initial 
solution estimate y n +i at the end of one time step, 

Vn+i = y n + /(t„,2/„)At + f s (r n ,y n ) ■ (asw n ). (23) 

The effect of the random fields is contained in the last 
term. The factor a s w n replaces the time- integral of the 
stochastic magnetic fields. For each site I of the array, 
there is a triple of unit variance, zero mean random num- 
bers (wf. 



ln> W ln' w l 



J produced by a random number gen- 
erator. The physical variance a s needed in the stochastic 
fields is defined by an equilibrium average over the time 
step. For an individual component at one site, that is 




dr' (K(r)K(r')) (24) 



When the FD theorem is applied to one dipole, this gives 
the variance of the random fields, 



a s = <t s (At) = V2aT At. 



(25) 



Thus, the individual stochastic field components, inte- 
grated over the time step, are replaced by random num- 
bers of zero mean with the variance o~ s ■ 

Then, in the corrector stage, the points y n and y n +i are 
used to get better estimates of the slope of the solution. 
Their average effect is used in the trapezoid correction 
stage: 

Vn+l = Vn + l; [f{r n ,Vn) + f(r n +i,Vn+i)] At (26) 

+ t; [fs(r n ,y n ) + fs(r n +i,Vn+i)] ■ (v s w n ). 

The error varies as C((Ar) 3 ), hence it is a second order 
scheme. It is important to note that the same random 
numbers w n are used in this corrector stage as those ap- 
plied in the predictor stage, for this individual time step. 

In the actual implementation of this method, the func- 
tions / and f s do not need to be specifically identified. 



Instead, it is seen that the change in any spin over a time 
step, A/2 = J dr-^jl, depends linearly on hAr from the 
deterministic part of the differential equation, and lin- 
early on J dr h s (r) from the stochastic part. But the 
stochastic contribution is replaced by random numbers 
of the correct variance, that is, 



hAr 



(27) 



Thus, the Euler predictor step is efficiently carried out 
by evaluating the combined deterministic plus stochastic 
field contributions, for an individual site, like 



fi = (i + Ajx, 
A/2 = A x W- a (A x 9)] 



(28) 
(29) 



where the needed effective field that updates this site is 
the deterministic/stochastic combination, 



g = h[fi] At + a s w 



(30) 



This latter relation generates an effective field, based 
on the deterministic local magnetic field from expres- 
sion (|I2[) . combined with a stochastic fluctuating impulse. 
The same type of combination applies in the trapezoid 
corrector step. The updating field at the end of the time 
step is calculated, using the predicted position, together 
with the same random field, 

g = h[p] At + a s w (31) 
That leads to a different estimate for the spin change, 



A/t = /2x g-a(fixg) 



(32) 



Then the corrector stage gives the updated spin accord- 
ing to their average 



A(t + At) = /2(r) 



A/2 



(33) 



This algorithm does not ensure the conservation of spin 
length. Thus, the length of /2 can be rescaled to unity 
after the step. However, in practice, when that be- 
comes necessary, the scheme is already becoming unsta- 
ble, which can happen if the net field \g\ is greater than 
1. The best way to avoid numerical instability is to make 
sure the time step is sufficiently small, so that \g\ <C 1, 
or even better, \Ajl\ <C 1, to the precision desired. 

The integration requires a sequence of quasi-random 
numbers (the w n stochastic fields ) with a long period, 
so that the simulation time does not surpass the pe- 
riod of the random numbers. We have used the gnerator 
mzranl3 due to Marsaglia and Zaman 2 - , implemented in 
the C-language for long integers. This generator is very 
simple and fast and has a period of about 2 125 , due to 
the combination of two separate generators with periods 
of 2 32 and 2 95 . 
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B. The dipole fields 

The calculation of the dipole term in the local magnetic 
field, (fT2")l . easily consumes most of the calculational ef- 
fort. In order to have correct results, the fields at all 
dipole sites must be calculated both in the predictor and 
corrector stages of a step. We consider a system with 
open boundaries, rather than periodic boundary condi- 
tions. In this way, the dipole calculation is simpler; there 
are N(N — l)/2 dipole field contributions to be found at 
any time. Some gain in speed can be achieved by com- 
puting tables of the pair displacements f{j and the powers 
rfj that are being re-used during the calculations. 

One of the best ways to speed up the calculation of 
the dipole fields for larger systems is to write their cal- 
culation as a convolution of a Green's function with the 
source dipoles, and calculate that convolution in recipro- 
cal space, transforming between real and reciprocal space 
with a fast Fourier transform (FFT). This is simplest if 
the locations of the dipoles form a simple square grid. 
Although that is true here, that grid is rotated at 45° 
to the x-axis, which brings some complications and only 
works for the square lattice spin ice. 

Instead, we consider that the spin ice involves unit cells 
on a square lattice, where each cell has a two-atom basis 
(i.e., two-dipole basis). For the cell whose lower left cor- 
ner is at position r k = (x k ,y k ) = ( , mk,nk)a, define the 



two dipoles present. On the "1" sublattice, 

//i(fk) at ffci = (x fc + \a,y k ), (34) 

and on the "2" sublattice, 

lh{rk) at f%2 = (xk,Vk + 50). (35) 

If there is an arbitrary source dipole fl at the origin, then 
the dipolar field h d it creates at position r — (x,y,z), 
according to the first term in (fl2|) . is well-known, 

, (2x 2 -y 2 ixy \ /V\ 
h d (f) = -,[ 3xy 2y 2 - x 2 • \ ^ \ (36) 
r \ -r 2 j V M * j 

(with all distances measured in lattice constants). This 
can be used to get the field produced from either sublat- 
tice. Summing over source dipoles, the 3x3 matrix is a 
Green's operator G(r) acting on the dipoles at discrete 
lattice sites. (Here the tilde is used only to indicate a 
3x3 matrix quantity.) However, to account for the two- 
atom basis, the Green's matrix is expanded to have an 
extra pair of indices that refer to the sublattice, one for 
the field point (a) and one for the source point (/J). The 
Green's matrix for the field produced at point r ka due to 
the source dipole at point rip is 



G a p(r k ,n) 



\r ko 



2(x ka - xip) 2 - (y ka - yip) 2 3(x ka - xip)(y ka - yip) 
3(x ka - xip)(y ka - yip) 2(y ka - yip) 2 - (x ka - x t p) 2 







\r ka - rip\ 



(37) 



It will be important to keep in mind that G a p actually 
depends only on the differences of the unit cell positions, 
r k i = f k — fi- Now the dipole field on the a sublattice, 
for the cell at r k , is given by a discrete convolution 



tices (a — l,/3 = 2) has a Cartesian element, 



N c 2 

1=1 p=\ 



Up{n) 



(38) 



Using fairly obvious notation, pp(fi) is the dipole on the 
j3 sublattice for the unit cell at rj. The dot operation 
represents the matrix multiplication, i.e., an implicit sum 
over the Cartesian components of G a p and pip. Written 
this way, the same formula could apply to other lattices 
of interest, such as honeycomb, Kagome, etc. Note that 
even at r k — f\ =0, there are contributions that must be 
included, corresponding to the interactions between the 
sublattices within an individual unit cell. For a specific 
example using the square ice sites, one can see that one 
particular interaction involving the two different sublat- 



G^(f k ,fi) = d- 



2(x k -xi + %) 2 



(Vk -vi- §) 2 



[(^.-X i + f)2 + (y fc - 2/i -f)2] 5 / 2 

This is nonzero when f k — rj. Also, the element 
G 21 (r k ,ri) with source and observer sublattices inter- 
changed can be obtained by changing the sign of a; they 
are not the same. A similar term with source and ob- 
server on the same sublattice is 

G n( r fe,n) = a gyj (40) 

{(x k - xi) 2 + (y k - yi) 2 } 

This is equal to G^{r k ,ri). It is divergent at r k = rj, 
however, that would be an interaction of a dipole with 
itself, a term that must be excluded by definition. 

G a p depends only on the displacements, f k i = r k — fi, 
which form another square lattice. Then one can find its 
Fourier transform, using a fast Fourier transform (FFT), 
setting the arbitrary source point to the origin. The 
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Fourier transform of jlp is also determined. The convolu- 
tion in real space becomes a simple product of G a p and 
ftp in Fourier space, which can then be transformed back 

to real space by an inverse FFT to obtain h d . Although 
there is considerable overhead, for larger systems the 
speedup is tremendous (N In N operations) when com- 
pared to doing the N sums with N terms to get the local 
dipole fields. 

To apply the simplest FFT method, the size of the 
grid of primitive cells must be 2 x x 2 » with integers 
N x and N y . To avoid the wraparound problem, so that 
the system being simulated is really a single copy of the 
desired L x x L y size, one needs to choose N x and N y 
large enough so that 2 W - > 2L X , and 2 N « > 2L y . This 
ensures that the periodic copies of the system, inherent 
in the application of the Fourier transform, do not "see" 
or interfere with each other in the convolution. 

It can also be kept in mind that there are some sym- 
metries that reduce the calculational overhead. Displace- 
ments only on the 1-sublattice or only on the 2-sublattice 
are the same, so for any of its Cartesian components, 

Gn = G 22 . (41) 

Also, the matrix is symmetric in the Cartesian indices, 
for any sublattice indices, 

G% = G%. (42) 

Furthermore, the interactions with both source and ob- 
server on the same sublattice are symmetrical in their 
interchange, 

G n (fki) = G n (rik) = G 22 {f kl ) = G 22 (r ;fc ) . (43) 

Then this implies their Fourier transforms are pure 
real, leading to some reduction in the computations 
needed. On the other hand, there is no symmetry gener- 
ally between different sublattices on different unit cells, 
so Gi 2 (rfc/) 7^ G 2 i(r*fc/), see the discussion after Eq. 
(|3"9")l . But the interactions within a cell must avoid self- 
interactions, so we do define Gn(0) = G 22 (0) = 0. The 
interactions between sublattices on the same cell depend 
only on squared displacements, so Gi 2 (0) = G 2 i(0) ^ 0. 
For fki 0, these symmetries result in 12 independent 
elements in G(fki) (each of Gn, Gi 2 , and G 2 i have 4 
independent elements), in contrast to the 4 independent 
elements needed for a single sublattice, see the matrix in 
Eq. (23). 

IV. THERMAL EQUILIBRIUM PROPERTIES 

Mostly the magnetic properties of spin-ice materials 
are investigated in an approximation of zero tempera- 
ture, because the fundamental interaction strengths of 
the anisotropy energies and the dipolar energies are much 
greater than UbT at room temperature (Models A,B,C). 



Even so, there could be energetic thermal fluctuations 
in a magnetic system even at low temperatures, in any 
situation where the magnetic fluctuations are large, such 
as near a reversal point in a hysteresis loop. This might 
possibly lead to enhancement of specific heat in such a 
situation, and of course, thermal rounding of the reversal 
paths in magnetization hysteresis loops. Thus it could 
be interesting to have some calculations of the energy, 
specific heat, and also of magnetic susceptibilities in a 
situation of thermal equilibrium. 

The time evolution from the Langevin dynamics can be 
used to get thermal averages, as an alternative to Monte 
Carlo calculations. As long as the simulation time is long 
compared to any physical relaxation time, a sequence of 
energy samples E n and total system magnetization sam- 
ples M n can be averaged, and their fluctuations can be 
used to estimate the specific heat and magnetic suscep- 
tibility. Suppose there are N s samples taken from the 
time evolution. The average energy for all N islands is 
estimated as 

( E ) = W^ En (44) 

s n 

with a "measurement error" estimated from its standard 
deviation and the number of samples, 

A£=-^|=, of = (E 2 ) (E) 2 . (45) 

The total heat capacity of the system is determined from 
the fluctuations in the energy, 

C N = /3 2 {(E-{E)f) = p 2 a E , (46) 

where /? = (fc^T) -1 , and from that we obtain the specific 
heat per island, c = Cn/N. The error in the heat capac- 
ity really needs to be estimated to be sure that the "mea- 
surement" is reliable. It can be calculated by finding the 
standard deviation of the quantity z = (E — (E)) 2 from 
which Cm was obtained, which turns out to be found 
from the following averages, 

a\ = (E 4 ) - (E 2 ) 2 +4(E) [2(E 2 )(E) - (E 3 ) - (E) 3 ] . 

Then the error in Cjv is 

AC = /3 2 4^ (48) 

and the error in specific heat per island is Ac = AC/N. 
Similar relations were used for the thermally averaged to- 
tal magnetic moment of the system and the components 
of magnetic susceptibility. For instance, the susceptibil- 
ity per island, Xxxi is found from fluctuations in total 
magnetic moment of the system, M x = /if , 

Xxx - ^ ((M x - (M x )f) = ^a 2 Mi , (49) 

and its error A\ xx is found from relations similar to those 
in (|47p and (l4"81 used for the heat capacity error, but with 
E replaced by M x . 
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Due to the fluctuations caused by the temperature in 
the simulations, the calculations of C and \ arc generally 
not as precise as those of (E) and (M) , without making 
very long runs. Especially as mentioned above, these 
calculations are difficult in any physical situation where 
the magnetization is on the verge of reversal, where the 
fluctuations are greatest. 

For these thermodynamics simulations, the system was 
started in a random state, with the temperature ini- 
tially set at the highest value in the range of interest. 
For the chosen temperature, data samples were taken 
at some appropriate time interval, that depends some- 
what on the energy couplings. For coupling parameters 
ki,k% on the order of 0.2 or less, and d several orders 
smaller (Models A,B,C), we found that a Heun time step 
At = 0.01 was sufficient to insure proper energy con- 
servation at zero temperature. Using this time step for 
finite temperature together with a — 0.1, we averaged 
over N s — 4000 data samples separated by sampling time 
interval Ar s = 10 3 At = 10.0 . An initial time interval 
corresponding to 100 samples was allowed for relaxation 
before samples were taken. Simulations would be left to 
run even longer than 4000 samples, if necessary, until the 
percent error in the magnitude of the total system mag- 
netization was found to be less than 0.1%. The final state 
at one temperature was then used as the initial state for 
the next lower temperature in the calculation. 

For models D and E, the dipolar coupling is much 
stronger, and this requires a smaller Heun time step, 
At = 0.001, to insure proper dynamics at finite tem- 
perature and energy conservation at zero temperature. 
Besides this, the calculation parameters for averages were 
the same as for models A,B,C, e.g., N s = 4000 and taking 
samples at sampling time interval Ar s = 10 3 Ar, while 
waiting for 0.1% or better precision in the system mag- 
netization. 

Some thermodynamic results for the Wang et al. parti- 
cles (Model A) are shown in Fig. [2j versus scaled temper- 
ature T ■ A 16 x 16 grid of cells was used (N = 2 x 16 2 = 
512). As the dimensionless energy coupling constants 
are small numbers, the only interesting effects are ob- 
served for T < 0.1 . Near T ~ 0.02 there are peaks 
in specific heat and in the in-plane components of mag- 
netic susceptibility. As mentioned earlier, though, these 
features would appear only greatly above room tempera- 
ture, which is marked with arrows. At these "high" tem- 
peratures, one would expect other modifications would 
take place (besides magnetic effects) and the model would 
not be applicable. Note that the monopole charge den- 
sity in Fig. does not go to zero at very low tempera- 
ture here. It does, however, make a transition to a lower 
value. This is an indication that the system did not find 
a state close to a ground state. There is frozen-in disor- 
der at lower temperatures. It is also an indication that 
the time scale for thermal relaxation to an equilibrium 
configuration was longer than the time interval used for 
averaging. Both the discrete and continuous definitions 
of p m exhibit similar behaviors. The order parameter 
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FIG. 2: (Model A) For a 16 x 16 grid of particles as used in 
Wang et al. with indicated parameters, (a) the internal energy 
and specific heat per site versus scaled temperature; (b) the 
components of the magnetic susceptibility at zero external 
field; (c) The monopole density, Eq. ([8} as determined from 
discrete and continuous charge definitions, Eqs. ([6| and 0. 
The vertical arrows show room temperature: it is essentially 
unaccessible in this dynamics. 
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Model B16: [1=100 nm x 12.5 nm x 5 run, a=200 nm - 
d=4.87e-5, k, =0.146, k =0.164, £=4.56 aJ - 
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Model CI 6: [1=40 nm x 80 nm x 4 nm, a=80 nm 
d=1.56e-4, k, =0.148, k =0.120, £=0.934 aJ 
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FIG. 3: (Model B) For a 16 x 16 grid of smaller sized particles 
with indicated parameters and a lower energy scale, (a) the 
average internal energy per site and the specific heat per site 
versus scaled temperature; (b) the components of the mag- 
netic susceptibility at zero external field; (c) the monopole 
charge density. 



FIG. 4: (Model C) For a 16 x 16 grid of still smaller sized 
particles with indicated parameters and an even lower energy 
scale, (a) the average internal energy per site and the specific 
heat per site versus scaled temperature; (b) the components 
of the magnetic susceptibility at zero external field; (c) the 
monopole charge density. 
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Z (not shown) stayed close to zero for the whole tem- 
perature range shown. That is further indication of the 
system staying far from a ground state, where it would 
have reached one of the values ±1. 

Similar results for Models B and C, for smaller parti- 
cles, are shown in Figures [3] and 0] These suggest that for 
the typical square lattice spin-ice using Py as the mate- 
rial, the room-temperature thermodynamics is nearly the 
same as that at zero temperature. There would be some 
limiting specific heat C « ks and non-zero value for the 
in-plane susceptibility. However, the monopole density 
has extreme difficulty to go to zero while scanning from 
high to low temperature. Then, in fact, the dynamics is 
a low temperature dynamics in a disordered non-ground 
state (and non-equilibrium) configuration. 

The thermodynamic results for 16 x 16 theoretical 
Models D and E are shown in Figures [5] and [6] Their gen- 
eral structure is similar: Model D has a strong peak in 
specific heat near T ~ 0.22 and a more rounded peak in 
Xxx ~ Xyy at a slightly higher temperature. The peak in 
specific heat in Model E is closer to T ~ 0.23, and again 
the susceptibility peak is centered a little higher in tem- 
perature. For all of the models studied, the out-of-plane 
magnetic susceptibility Xzz is considerably smaller than 
Xxx, and there is only a weak temperature dependence. 
Models D and E, notably, do reach thermodynamic equi- 
librium in the simulations. This is seen clearly in the 
plots of the order parameters Z and p m . Now p m , for 
both discrete and continuous forms, tends towards zero at 
low temperature, as expected for the system moving to- 
wards a ground state. Further, the ground state overlap 
order parameter, Z, tends to go towards unity as T — > 0; 
this is the strongest indication of approaching one of the 
ground states. At higher temperatures above the peaks 
in C and x we see that Z becomes quite close to zero; the 
system is more random and far from a ground state. In 
the same high-temperature region, the monopole density 
tends towards some positive constant limiting value. The 
discrete definition for p m apparently undergoes a stronger 
change in value as the system makes a transition from its 
low to high temperature behavior. 



V. HYSTERESIS CALCULATIONS 

A simple experiment to investigate the magnetic prop- 
erties of spin ice is the response in an applied external 
field (hysteresis calculation) . To get a general impression 
of the physical response in any square spin ice, hysteresis 
calculations were carried out for Model D at fixed scaled 
temperature T = 0.1 . These were calculated the same 
as for the thermodynamics, except that it was adequate 
to average over shorter sequences, N s = 1000, at each 
step of the applied field. The system was initially set in 
a random configuration, but with the maximum positive 
applied field. The field was scanned to lower and nega- 
tive values along some axis (either x or at 45° to +x) and 
then allowed to come back to the starting value. In order 
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FIG. 5: (Model D) For the model with D = Ki = K 3 = ^e, 
(a) the average internal energy per site (in energy units e) 
and the resulting specific heat per site versus scaled temper- 
ature; (b) the components of the magnetic susceptibility at 
zero external field; (c) the monopole charge density together 
with the ground state overlap order parameter Z. 
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FIG. 6: (Model E) For the alternative model with D = Ki = 
YqC, K3 = ie, (a) the average internal energy per site (in 
energy units e) and the resulting specific heat per site ver- 
sus scaled temperature; (b) the components of the magnetic 
susceptibility at zero external field; (c) the monopole charge 
density together with the ground state overlap order param- 
eter Z. 



FIG. 7: For the model with D = Li = K3 — -^e, at tem- 
perature ksT = y-jE, (a) the averaged magnetization per site 
versus external magnetic field ft cxt applied along the x-axis; 
(b) the order parameter Z versus /i ex t- The field strength was 
initially set at /i, cxt = 0.8, and scanned to /i cxt = —0.8, then 
back to the starting value as in a hysteresis loop calculation. 

to interpret the results, it was also important to calculate 
the order parameter Z and the monopole charge density 
p m during the hysteresis scan. 

The results for field applied along the x axis are shown 
in Figures [7] and [3 In fact, at this temperature, the 
model does not exhibit any hysteresis: the magnetiza- 
tion per island is the same in backward and forward scans 
of /iext- However, the magnetization shows regions with 
distinctly different slopes. In Fig. [7p, one sees that the 
order parameter Z , however, tends to take on either val- 
ues close to zero, at strong applied field, or, values near 
Z ss ±1, at weaker field. We note that this temperature 
T = 0.1 is on the low side of the specific heat peak for 
zerro magnetic field. Then this shows that in the central 
region of the MH graph, the system falls into states that 
are close to the ground states. These states, however, 
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FIG. 8: For the model with D = Ki = K 3 = ^e, at tem- 
perature ksT = jqE, the averaged monopole density versus 
external magnetic field /i cxt applied along the x-axis, for the 
calculations in Fig. [7J The results of both the discrete and 
continuous charge definitions are compared here. 



are slightly modified due to tilting of some of the dipoles 
according to the field strength. Hence, there is close to a 
linear response with h ex t, as the dipoles on the 2 nd sub- 
lattice, which are nearly perpendicular to the field, get 
tilted by it. 

The variation of monopole charge density with applied 
field, for the simulation in Fig. [71 is shown in Fig. El Both 
the discrete and continuous definitions are displayed. 
The discrete definition has more dramatic changes (al- 
though, its errors, not shown, may be larger). Especially, 
p m tends to zero (or a small value for the continuous ver- 
sion) over the same applied field range where Z ±1. 
This confirms clearly that the central region of the MH 
graph corresponds to the system being in states close to 
the ground states. 

For applied field along an axis at 45° to the +i-axis, 
the situation is similar, see Figures [S] and [TO] The MH 
and ZH graphs are nearly the same as for those for ap- 
plied field along x. In this case, however, the applied 
field must be causing both sublattice dipoles to tilt at 
stronger fields. There is a difference, then, in the charge 
density plot, see Figure [TU] Again, in the central region 
near weak /i ex t, the monopole density tends to zero, as 
expected for the system being close to one of the ground 
states. At strong field, however, the discrete charge den- 
sity also tends towards zero (as does Z). In this case, 
though, as the dipoles all tend to align at 45° to ±x, 
the net number of dipoles pointing into any charge site is 
then forced to be zero. This clearly forces the monopole 
density found by the discrete definition to zero. One sees 
that the monopole density by the continuous definition, 
on the other hand, does not fall to zero at high field, and 
instead behaves the same as it does for field applied along 
x. 
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FIG. 9: For the model with D = Ki = K 3 = ie, at tem- 
perature ksT = y-jE, (a) the averaged magnetization per site 
versus external magnetic field /i ox t applied at 45° above the x- 
axis; (b) the order parameter Z versus /i oxt . The field strength 
was initially set at /i cxt = 0.8, and scanned to h cxt = —0.8, 
then back to the starting value as in a hysteresis loop calcu- 
lation. 



VI. DISCUSSION AND CONCLUSIONS 

We have studied the possibilities of spin dynamics in 
artificial spin ice systems. To do this, we have analyzed 
frustrated arrays consisting of two-dimensional square 
lattices of elongated magnetic nanoislands. The inter- 
nal structure of the magnetic nanoislands was taken into 
account in such a way that the dynamics of all spins in- 
side each magnet were inserted in the interactions. Then, 
depending on the island shapes, aspect ratios, sizes, el- 
ements and organization in the array, we have looked 
for possible departures from the Ising-like behavior com- 
monly found for these systems. In general, it is not easy 
to build square ice arrays with dynamics by using nanois- 
lands made of permalloy. As a consequence, we have 
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FIG. 10: For the model with D = ifi = Ka = j^e, at tem- 
perature ksT = jqE, the averaged monopole density versus 
external magnetic field h ex t applied at 45° above the rr-axis, 
for the calculations in Fig. [9] 

found that the systems without real dynamics (islands 
practically with an effective Ising behavior) have great 
difficulty to achieve the ground state (see, for instance, 
models A, B and C). Indeed, the order parameter Z (de- 
fined in Section II) never reaches the values —1 or +1 
(the two degenerate ground states), even for very low 
temperatures. This result agrees with all experimental 
studies&i^ concerning the square lattices. On the other 
hand, by considering fictitious materials (as described by 
the constants D, K\ and K 3 ), we found interesting de- 
viations from the Ising behavior and consequently, more 
easily thermalized spin dynamics along the array. For 
these systems (see models D and E), the ground state 
can be easily obtained for low temperatures. 

Some of the results obtained here can be directly com- 
pared to those of Ref. 12CJ, where the artificial square spin 
ice with point-like dipoles with Ising-like behavior was 
studied by using conventional Monte Carlo simulations. 
Here, for simplicity, the model which replaces the mag- 
netic moment of the islands by a point-like and Ising-like 
dipole (see, for instance, Refs.— ~ 16 i 20 ) will be referred to 
as model F. In model F, the ground state (Z — ±1) of 
the square ice can be easily obtained for very low tem- 
peratures, i.e., by using conventional Monte Carlo simu- 
lations; the ground state naturally appears at low enough 
temperatures. However, our results for models A, B and 
C, obtained by using a Langevin dynamics, showed that 
the ground state does not appear at low temperatures. 
This may indicate that there is a kind of dynamic con- 
straint (effectively, excessively long relaxation time) in 
the system's dynamics that prevents it from reaching the 
ground state over a moderate time of observation. In- 
deed, a similar result was found when the dynamics of 
the model F in the presence of external magnetic fields 
is considered 12 ' 27 " — . Thus, it may indicate that arti- 



ficial square spin ices made with Permalloy may never 
reach the ground state, since both external field dynam- 
ics and thermally driven dynamics have bottlenecks that 
prevent the access of the ground state. On the other 
hand, our results for models D and E, for materials with 
fictitious constants, indicate that it may be possible to 
access the ground state using another kind of material. 
It seems then that, by changing the island's anisotropies 
and interactions, the dynamical bottleneck is eliminated. 

This difference may be associated with the way that 
the system explores the phase space. First, let us con- 
sider systems where the islands have an Ising-like behav- 
ior. Since we are dealing with classical particles, each 
island must pass through an energy barrier to change its 
magnetization direction (even being Ising-like), either by 
the creation and propagation of a domain wall or by ro- 
tation of a single domain. Either way, there is no option 
for tunneling and some energy must flow to the island. 
This internal energy barrier was not taken into account in 
the Monte Carlo calculations of Refs i 14 ' 15 i 20 and probably 
this is why the ground state was obtained. The results 
of the present study and those from Refs.— 2Zr2& include 
the energy barrier for spin flips and we may expect that 
the existence of a huge energy barrier is responsible for 
the difficulty to access the ground state. Moreover, in 
models D and E, where the energy barrier is smaller, the 
ground state is accessible. It seems then that the key fac- 
tor blocking access to the ground state of artificial square 
spin ice is the energy barrier for spin flips. 

For model F, the specific heat exhibits a peak at a 
temperature around 7.2D/kB, where it is suspected that 
the string connecting the Nambu monopoles is broken 
and the system is able to support free monopoles 2 ^. The 
specific heat for models A, B, C, D and E also exhibits 
a characteristic peak. However, for models A, B and C, 
the peak arises for much higher temperatures (around 
20-D/fcs); it may seem unreasonable since the magnetic 
moment of the islands for these models has many more 
degrees of freedom than that of model F, Nonetheless, 
we remark that the ground state can not be obtained for 
models A, B, C even for zero temperature and, therefore, 
there is not a way of establishing equivalences between 
the results of these models and results of model F. Al- 
ternatively, the ground state of models D and E can be 
found at zero temperature and then, we can compare the 
results for the specific heat with the ones obtained for 
model F. Really, the peak in the specific heat for models 
D and E arises at a temperature around 2D/ks, which 
is about three times smaller than the case of model F. 
It is obviously expected because there are more spin de- 
grees of freedom for models D and E than for model F. 
Furthermore, the peak for model E occurs at a temper- 
ature slightly above that one for model D because the 
out-of-plane spin degree of freedom is more restricted in 
E. 

To complete this study, we have also calculated the 
hysteresis for the models that exhibit dynamics. It is an 
important calculation to get a general impression of the 
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physical response in any square spin ice, since the sim- 
ple experiments to investigate the magnetic properties of 
these systems is the response in an applied external field. 
Finally, we would like to say that we really hope that all 
the investigations developed here could help experimen- 
tal advances toward spin ice systems in which the ground 
state could be achieved and/or the transition rendered by 
appearance of free monopoles occurs around room tem- 
perature. Experimentally, a recent work was already ad- 
dressed in this direction. Indeed, Kapaklis et al~2- have 
proposed an experimental system (in an external mag- 
netic field) where thermal dynamics can be introduced by 
varying the temperature of the array. On a square lattice, 
they use a material (based on 5-doped Pd(Fe)) with an 
ordering temperature near room temperature to confirm 
a dynamical "pre-melting" of the artificial spin ice struc- 
ture at a temperature well below the intrinsic ordering 
temperature of the island material. Such a procedure is 



capable of creating a spin ice array that has real thermal 
dynamics of the artificial spins over an extended tem- 
perature ranged. The possibility of observing emergent 
monopoles is therefore conceivable, following the general 
approach that the authors of Ref. [l9| described in the 
design of spin ice arrays. This is a first step towards re- 
alization of artificial spin ices as conceived in models D 
and E. 
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